function [TSFC,T,Wf_dot] = engine_model(ac,alt,V_tas,throttle,ISA)

%SI units

% % % % % % %%% Test Only 
% % % engine_load(ISA);
% % % % % %%%

th=throttle;   %Aircraft throttle form idle [%]
speed=V_tas; %Aircraft speed [m/s]

flag=1; %-> temp means delta ISA in K
temp=ISA;

%Atmospheric properties, SI units

[p,rho,mu,T,sound,delta_isa]=isa_cos(alt,temp,flag);

%Calculate Mach number
M=speed/(sound);

%Load Engine parameters saved in engine.mat file
load('engine');
%% Interpolate engine database with current aircraft conditions

T = interp3(M_range,alt_range,th_range,T_matrix,M,alt,th); %Engine thrust in N;
TSFC = interp3(M_range,alt_range,th_range,TSFC_matrix,M,alt,th)/1000; %Engine TSFC in kg/N.s;
Wf_dot = (T*TSFC); %Engine Fuel flow in kg/s
